中国-东盟1 km分辨率NDVI数据集（2013）

李  静，曾也鲁，柳钦火*，仲  波，吴善龙，彭菁菁

中国科学院遥感与数字地球研究所遥感科学国家重点实验室，北京 100101

摘  要：中国-东盟1 km分辨率NDVI数据集（2013）是在5种1 km分辨率卫星遥感数据（Terra/MODIS、Aqua/MODIS、NOAA18/AVHRR、FY3A/VIRR和FY3B/VIRR）基础上，参照标准的MODIS-NDVI合成算法生产的、5-day分辨率的数据集。本数据集的合成算法分两步：首先，在5-day合成期内，利用基于带NDVI权重的核系数稳健拟合法，把多源数据按照残差阈值10%、20%和大于20%划分为三个质量等级，其中残差阈值大于20%的三级数据质量最差，不参与合成。然后，根据5-day合成期内各质量等级的有效观测数目，分别采取主算法或备用算法进行NDVI合成，获得了2013年中国-东盟时间分辨率为5-day、空间分辨率为1 km的NDVI时间序列。在中国黑河流域利用高分辨率ASTER影像，将本数据集与标准MODIS NDVI产品进行对比分析。结果表明，基于多源数据集生成的NDVI与ASTER-NDVI的决定系数R2显著高于基于单源数据集生成的标准MODIS-NDVI。该数据集采用Sinusoidal Tile Grid分幅，共40景，每景覆盖的经纬度范围为10°×10°。该数据集存储为.tif 格式，数据量为31.17 GB，被压缩为40个.zip文件，压缩后的数据量为9.89 GB。

关键词：NDVI；中国-东盟；二向性反射分布函数（BRDF）；核驱动模型；质量分级

DOI: 10.3974/geodp.2017.03.03

 

1  前言

植被指数是根据植被的光谱反射与吸收特性，将卫星遥感数据可见光与近红外波段的反射率进行线性或非线性组合，形成对地表植被变化敏感、可以有效反映植被覆盖程度的一系列指示因子的总称。植被指数与叶面积指数、净初级生产力等生物物理参数有着密切的关联，是植被覆盖的动态监测、成图和资源管理中的一个重要参数。归一化差分植被指数（Normalized Difference Vegetation Index, NDVI）是迄今为止应用最为广泛的一种植被指数。现有的全球尺度NDVI植被指数产品均是基于单一卫星传感器。由于受云、气溶胶等天气状况的限制，单一传感器在合成周期内提供的有效观测数目有限，这导致NDVI产品的合成周期过长。作为描述植被特征的植被指数，过长的合成周期，一方面会导致NDVI植被指数合成精度变低；另一方面，无法有效准确的识别出植被动态变化，这限制了NDVI合成植被指数产品对植被的指示作用。确定合成周期，最主要的考虑因素是保证一定的有效观测个数以确保合成精度。目前，基于单一传感器的合成方式，合成周期最优为10-day。充分利用目前在轨的多卫星观测，提高单位时间内的观测次数，是提高合成周期的有效途径。

中国-东盟1 km分辨率NDVI数据集（2013）[1]（MuSyQ-NDVI-1km-2013）由中华人民共和国科学技术部项目资助完成。本产品是基于目前可以方便获取的1 km分辨率全球覆盖的传感器数据产品研发的，具体包括Terra/MODIS、Aqua/MODIS、NOAA18/AVHRR、FY3A/VIRR、FY3B/VIRR等传感器数据。在对多源数据集、现有的合成植被指数产品以及算法进行分析的基础上，提出了基于多源数据集研发1 km分辨率5-day周期的全球合成植被指数产品算法体系。多源数据集可以在有限时间内提供比单一传感器更多的角度和更多次的观测。但是，由于传感器在轨运行时间及性能差异，多源数据集的观测质量参差不齐。因此，为更有效地利用多源数据集，算法体系首先对多源数据集进行了质量分级，根据观测合理性分为一级数据、二级数据、三级数据。三级数据为受薄云污染的观测，不用于计算。植被指数合成算法基本采用MODIS的植被指数合成算法，即基于半经验Walthall模型的BRDF角度归一化方法、CV-MVC法和MVC法的算法体系。利用该算法体系，分别对一级数据、二级数据计算合成植被指数，并进行质量标识。本产品在中国黑河流域与高分辨率卫星影像及目前主流的MODIS NDVI产品进行了对比分析。

2  数据集元数据简介

中国-东盟1 km分辨率NDVI数据集（2013）的名称、作者、地理区域、数据年代、时间分辨率、空间分辨率、数据集组成、数据出版与共享服务平台、数据共享政策等信息见表1。

3  数据研发方法

3.1  算法原理

产品算法的目标是充分利用多源遥感数据集提供的多角度多时相的优质观测，提高合成植被指数的时间分辨率和精度。由于传感器的观测噪声、残云存在等因素的影响，预处理后的多角度数据集会存在一些噪声及误差较大的观测值。因此，本算法设计了基于观测质量分级的植被指数合成算法体系。根据质量分级之后各个级别观测数目的多少，分别采取主算法与备用算法进行多源遥感数据集NDVI植被指数合成。

3.1.1  多源遥感数据质量分级

如图1所示，首先对多源数据组成的多角度输入数据进行质量分级。受残云或大气的影响NDVI值会偏小，因此NDVI可以作为检验反射率产品质量的指示因子之一。因为NDVI随太阳入射、观测角度变化显著，所以在对不同入射/观测角度下的NDVI进行比较时，引入核驱动模型消除NDVI角度效应，归一化到天顶方向观测的等效NDVI之后进行比较。多源遥感数据集质量分级包括以下4个步骤：

步骤一：选定一性能最佳的传感器(如MODIS)作为基准，对多源数据集中来自不同传感器的反射率数据进行波段转换，将所有反射率数据投影到相同的坐标系下，进行云检测等质量检查。

表1  中国-东盟1 km分辨率NDVI数据集（2013）元数据简表

条  目

描  述

数据集名称

中国-东盟1 km分辨率NDVI数据集（2013）

数据集短名

MuSyQ-NDVI-1km-2013

作者信息

李静 R-7298-2016, 中国科学院遥感与数字地球研究所遥感科学国家重点实验室, lijing01@radi.ac.cn

曾也鲁G-6465-2017, 中国科学院遥感与数字地球研究所遥感科学国家重点实验室, zengyl@radi.ac.cn

柳钦火S-1647-2016, 中国科学院遥感与数字地球研究所遥感科学国家重点实验室, liuqh@radi.ac.cn

仲波 L-4528-2016, 中国科学院遥感与数字地球研究所遥感科学国家重点实验室, zhongbo@radi.ac.cn

吴善龙S-1572-2016, 中国科学院遥感与数字地球研究所遥感科学国家重点实验室,wsl0579@163.com

彭菁菁S-1261-2016,中国科学院遥感与数字地球研究所遥感科学国家重点实验室, rspeggy@163.com

地理区域

地理范围：10°55′12″S-53°32′24″N，73°37′12″E-141°0′36″E。包括：中国、印度尼西亚、马来西亚、菲律宾、新加坡、泰国、文莱、越南、老挝、缅甸和柬埔寨

数据年代

2013年

空间分辨率

1 km

时间分辨率

5 day

数据格式

.tif

数据量

9.89 GB（压缩后）

数据集组成

数据集由40个压缩文件组成。每个文件解压后为该景2013年73个时相的数据，每个时相包括4个文件：NDVI.tif、NDVIqc.tif、EVI.tif和EVIqc.tif，分别代表：NDVI、NDVI质量标识、EVI和EVI质量标识

基金项目

中华人民共和国科学技术部 ( 2012AA12A304, 2013AA12A301)

出版与共享服务平台

全球变化科学研究数据出版系统 http://www.geodoi.ac.cn

地址

北京市朝阳区大屯路甲11号100101，中国科学院地理科学与资源研究所

数据共享政策

全球变化科学研究数据出版系统的“数据”包括元数据（中英文）、实体数据（中英文）和通过《全球变化数据学报》（中英文）发表的数据论文。其共享政策如下：（1）“数据”以最便利的方式通过互联网系统免费向全社会开放，用户免费浏览、免费下载；（2）最终用户使用“数据”需要按照引用格式在参考文献或适当的位置标注数据来源；（3）增值服务用户或以任何形式散发和传播（包括通过计算机服务器）“数据”的用户需要与《全球变化数据学报》（中英文）编辑部签署书面协议，获得许可；（4）摘取“数据”中的部分记录创作新数据的作者需要遵循10%引用原则，即从本数据集中摘取的数据记录少于新数据集总记录量的10%，同时需要对摘取的数据记录标注数据来源[2]

 

步骤二：若经过质量检查的总观测数目N不少于5个，则根据各观测计算的NDVI均值减去0.3作为阈值，小于阈值是明显异常的观测值，直接作为三级数据；若确定三级数据之后的其余观测数目N仍不少于5个，则进入步骤三；若N小于5个，则进入步骤四。

步骤三：若经过步骤二剔除异常观测值后的MODIS有效观测数目M不少于5个，则将其作为数据输入，采用基于稳健估计的核驱动模型拟合出核系数；若经过步骤二剔除异常观测值后的MODIS有效观测数目M少于5个，则直接采用对应像元、对应时相的MODIS BRDF产品，获取相应的核系数；通过核系数对各观测量方向的NDVI进行角度校正，得到天顶方向的等效观测NDVI；选取等效观测NDVI的次大值作为基准，以相对误差的10%、 20%为界划分出一、二、三级数据。

图 1  多源遥感数据集质量分级流程图

步骤四：若经过步骤一质量检查的总观测数目N少于5个，或经过步骤二剔除异常观测值后的总观测数目N小于5个，则直接采用对应像元、对应时相的MODIS BRDF产品，获取相应的核系数；通过核系数对各观测量方向的NDVI进行角度校正，得到该观测在天顶方向的等效观测NDVI，选取等效观测NDVI的次大值作为基准，以相对误差的20%为界划分出二、三级数据。

3.1.2 主算法——基于BRDF模型的拟合法

随着噪声移除技术尤其是BRDF校正技术发展，植被指数角度归一化合成方法得以发展。它采用BRDF模型将合成周期内所有无云观测值逐波段逐像元拟合至某一照射条件下星下点等效反射率值，再计算得到合成值[3–5]。变化的入射/观测几何是植被指数方向效应的主要来源，且生物物理参数如LAI、FPAR等的反演模型都是基于星下点植被指数值。角度归一化合成方法基于辐射传输模型反演，移除变化的太阳-目标-传感器几何条件影响，代表了合成算法的主要进步和发展方向。Leeuwen等[6]对比了多种植被指数合成方式，发现Walthall BRDF模型对于许多植被类型计算星下点等效植被指数效果最好。因此，当一级数据的观测数目不少于5个，或一级数据与二级数据的观测数目之和不少于5个时，采用Walthall模型将反射率标准化到星下点并计算以星下点为基准的植被指数[7]，BRDF模式订正公式如式（1）所示：

             (1)

式中，ρλ是大气订正的反射率；θv是卫星天顶角；φs是太阳方位角；φv是卫星方位角。模型参数aλ，bλ和cλ用最小二乘法拟合得到。NDVI植被指数合成公式如式（2）所示：

                         (2)

3.1.3  备用算法

（1）平均合成MC法：MC（Mean Compositing）法计算合成周期内各观测的平均值，基于平均反射率值计算的植被指数作为合成植被指数[8]。当一级数据与二级数据的观测数目之和为2-4个时采取该算法。采用MC平均合成法，既有利于进一步减少传感器之间的差异，又能够充分利用有效观测信息避免有效信息的浪费，同时取平均的做法使得合成结果更加稳定。

（2）直接计算VI法：采取观测的反射率值直接计算植被指数作为合成植被指数。当一级数据与二级数据的观测数目之和为1个时采取该算法。该算法对地表非朗伯特性和观测几何变化考虑不足，易受残云、薄雾、传感器噪声等因素影响，稳定性稍差。

（3） 最大合成MVC法：MVC（Maximum Value Composite）法以合成周期内所有观测的最大植被指数作为合成植被指数。当不含一级数据与二级数据时采取该算法。MVC法在设计的时候重点考虑大气的影响，合成值更接近星下点观测值，且运算效率高，不足之处在于没有考虑地表非朗伯特性。

3.2  技术路线

如图2所示，在本产品算法中，基于不同的观测数据个数及数据质量情况，使用了不同的合成算法。算法边界的判断依据为合成周期内经过预处理后的无云观测的观测数据个数及质量。假设预处理后的无云观测数为N，经过数据质量分级后，一级质量数据为L1，二级质量数据为L2，三级质量数据为L3，N=L1+L2+L3。算法具体流程为：

（1） 当L1≥5时，用一级质量数据基于Walthall模型拟合法合成植被指数。此条件下，多传感器提供了质量较高、一致性较好的多时相观测，并采用合成精度较高的BRDF模型拟合法，是具有最优精度的合成结果。质量描述QA设为0。

（2） 当L1<5且L1+L2≥5时，用一级和二级质量数据基于Walthall模型拟合法合成植被指数。此条件下，具有较好一致性的一级数据个数不足以拟合Walthall模型的系数，将有一定观测误差存在的二级数据和一级数据一起应用于Walthall模型拟合法。由于观测数据的不一致性增大，合成精度会降低。质量描述QA设为1。

（3） 当1<L1+L2<5时，由于多角度观测数据不足，无法应用Walthall模型拟合法，采用备用算法的MC法计算合成植被指数。MC法先计算各观测的反射率的平均值，再计算植被指数。此条件下，合成结果仍具有一定的角度效应，其大小取决于观测角的分布，具有一定随机性。质量描述设为2。

图2  产品研发的技术路线

 

图3  中国-东盟1 km/5day分辨率NDVI植被
指数产品（审图号：GS(2015)1527号）

（4） 当L1+L2=1时，即一级二级质量的观测只有1个，此时采用VI法。质量描述设为3。

（5） 当L1+L2=0且N>0时，即没有一级二级质量的观测，只有三级质量观测时，采用MVC法。质量描述设为4。

（6） 当N=0时，即合成周期内没有无云观测，填充FillValue（FillValue = -999）。

4  数据结果与验证

4.1  数据集组成与空间分布

本数据集采用Sinusoidal Tile Grid分幅，共40景，每景覆盖范围为10°×10°。每景数据存放在一个文件夹里，包括该景2013年每5-day的共73个时相的数据。每个时相包括4个文件：NDVI.tif，NDVIqc.tif，EVI.tif，和EVIqc.tif，分别代表NDVI，NDVI质量标识EVI和EVI质量标识。

中国-东盟1 km分辨率NDVI植被指数合成结果如图3所示。由图3知，NDVI产品的空间连续性较好，数据缺失比例较低，稳定性较高，基本能满足对植被状态进行大范围监测的需求。

4.2  数据结果验证

4.2.1  直接验证

图4  中国黑河流域ASTER反射率数据假
彩色合成图 (DOY192)

如图4所示，将中国黑河流域分辨率为15 m的ASTER反射率数据经辐射定标，6S模型大气校正等预处理环节之后，经波段转换后投影到1 km分辨率NDVI产品对应的坐标系下，并聚合到1 km分辨率尺度，生成对应的植被指数，作为参考值与生产的产品进行直接验证。ASTER数据对应的地表覆盖类型主要为农田与荒漠，选取地块相对完整、成片连续的农田作为验证区域。验证区域大小为22 km×24 km，对应528个分辨率为1 km的像元。

图5分别列出了4种植被指数产品与ASTER数据生成的1 km分辨率NDVI参考图的对比结果图。4种植被指数产品分别为基于MODIS+FY多源数据集、基于MODIS单源数据集、基于FY单源数据集生产的产品，以及MOD13A2产品。决定系数（R2）可以表征生成的植被指数产品与植被指数参考图的相关性。由图5可知，相比于基于MODIS或FY单源数据，基于多源数据集生产的产品与NDVI参考图的R2分别从0.337或0.107提高到0.395，表明基于多源数据集生产的产品与NDVI参考图的相关性更好。其次，基于多源数据集生

图5  黑河流域农田站点NDVI产品与ASTER数据生成的参考图对比结果（DOY186）

产的产品均方根误差（RMSE）为0.102，拟合残差（Residual）为0.038，小于基于MODIS或FY单源数据生成的产品。

从整体上看，基于MODIS单源数据集生成的产品与NDVI参考图相关性较好，但拟合残差较大（0.042）；基于FY单源数据集生产的产品与NDVI参考图的RMSE较小，但拟合残差较大（0.044），且R2较低（0.107），与NDVI参考图的相关性较差。以上分析表明基于多源数据集生成的产品比基于单源数据集，精度较高，且更稳定。

与MOD13A2产品相比，基于多源数据集生成产品的RMSE和拟合残差与其相当，但R2高于MOD13A2产品。这表明基于多源数据集生成的产品与NDVI参考图的相关性更好。其次，MOD13A2产品的时间分辨率为16-day，而基于多源数据集生成产品的时间分辨率为5-day，显著提高了产品的时间分辨率。

4.2.2  交叉验证

交叉验证是一种间接的验证方式，是在没有地面测量数据的支持下，将不同产品统一到相同的投影坐标系和空间分辨率下，通过不同产品之间的相互比较来评估产品的效果。通过不同产品之间的交叉验证，评价产品的时空一致性，能在代表全球植被类型分布的大量站点进行。

（1）同一时相产品比较

在黑河流域选取地块相对完整、成片连续的农田作为验证区域。验证区域大小为22 km ×24 km，对应528个分辨率为1 km的像元。将生成5-day的NDVI产品与16-day的MOD13A2产品进行逐像元的对比分析，如图6所示。其中，基于多源数据集的NDVI产品生产时间区段为DOY193-197，DOY198-202，DOY203-207，一共3个生产时间区段；MOD13A2产品的对应生产时间区段为DOY193-208，一共1个生产时间区段。由图6可知，生产的NDVI植被指数产品与MOD13A2相比，均方根误差（RMSE）分别为0.048、0.072和0.055，拟合残差（Residual）分别为0.037、0.038和0.043，平均误差控制在0.1以内。基于多源数据集的植被指数合成产品将时间分辨率从MOD13A2的16-day提高到5-day。

（2）时间序列比较分析

合成植被指数产品的一个重要功能是监测植被的动态变化。在黑河流域农田区域随机选取5个像元，并与MOD13A2产品进行时间序列分析，如图7所示。其中，MOD13A2产品的起始时间分别为DOY177、193、209；基于多源数据集的产品起始时间分别为DOY183、188、193、198、203、208。由图7可知，从DOY177到DOY209，MOD13A2产品与基于多源数据集的合成植被指数产品都经历了先增加后减小的变化趋势。MOD13A2产品由于是在16-day时间长度上进行的合成，因而曲线的变化趋势显得更加平滑；基于多源数据集的产品因其是在5-day时间长度上进行的合成，因而曲线的变化幅度更加显著。从MOD13A2产品来看，5个站点的NDVI都在DOY193达到峰值。而从时间分辨率更高的多源NDVI合成产品来看，第3号站点与第5号站点在DOY193达到峰值，但第1、2、4号站点在DOY188就已经达到峰值。如果用MOD13A2的16-day合成产品，将难以捕捉到这一生长峰值的准确时间信息。表明更高的时间分辨率有利于提取植被的生长峰值信息与动态监测植被的长势变化。

图6  黑河流域农田站点基于多源数据集生成的NDVI产品与MOD13A2产品对比结果

图7  黑河流域5个农田站点基于多源数据集生成的
NDVI产品与MOD13A2产品时间序列对比结果

5  讨论和总结

随着目前全球在轨卫星的不断增多，同一天、甚至相近时相都有多传感器过境同一地区。同一天对同一地点的多个传感器观测组成了对地表植被不同太阳/观测平面的多角度数据集，比单一传感器具有更丰富的角度和观测信息，为提高植被指数合成产品的精度和时间分辨率提供了新的机会和可能。在此背景下，本算法设计了基于卫星组网构建的多源数据集生产全球1 km合成植被指数的产品算法。为了能够充分利用多传感器提供的优质观测提高植被指数合成精度及时间分辨率，本算法首先对多源数据集进行了质量分级。通过基于核驱动模型的迭代法，根据具有较好一致性的观测为优质观测的原则，将多源数据集分为三个质量等级。在质量分级基础上，对较优质量的观测，采用合成精度较高的基于BRDF模型拟合法；当数据质量变差、观测个数变少时，分别采用平均MC法、VI法以及MVC法，进而构建得到基于多源数据集的植被指数合成算法。

本算法在实施和产品生产过程中可能遇到的问题主要包括两点：一是输入参数的不确定性，即输入反射率数据的不确定性；二是多源数据集观测个数较少时，质量分级算法无法执行，经过异常值检验的所有观测都将参与植被指数合成，会为结果带来不确定性。

（1）输入参数的不确定性

多源遥感数据预处理过程存在不确定性，受到多源遥感数据交叉辐射定标精度、几何校正精度、大气校正精度、以及波段转换结果精度等的影响，造成反射率偏差，从而影响植被指数合成结果精度。云检测结果会影响合成算法的选择，对合成植被指数结果有很大影响。如果有云像元被检测为无云，尤其对于BRDF拟合，有云像元反射率影响植被特征波段的值，易造成合成结果的偏差。

（2）多源数据集观测个数较少时，合成结果有可能会出现异常

虽然算法体系中设计了质量检查过程，但是当观测数据较少时，质量分级算法将无法执行。数据质量检查仅仅为与最大NDVI的差值小于0.3。经过这个质量检查后的所有观测将直接参加植被指数合成的运算。当观测数据较少时，观测很有可能会受残云等影响而出现误差。由于无法做质量检查，合成结果有可能存在异常。

 

作者分工：李静、曾也鲁和柳钦火总体设计了数据集的研发流程；仲波、吴善龙和彭菁菁采集和处理了数据；李静和曾也鲁撰写了数据论文等。

参考文献

[1]       李静, 曾也鲁, 柳钦火等. 中国-东盟1 km分辨率NDVI数据集（2013）[DB/OL]. 全球变化科学研究数据出版系统, 2015. DOI: 10.3974/geodb.2015.01.16.V1.

[2]       全球变化科学研究数据出版系统. 全球变化科学研究数据共享政策[OL]. DOI: 10.3974/dp.policy.

2014.05 (2017年更新).

[3]       Duchenmin, B., Berthelot, B., Dedieu, G., et al. Normalisation of directional effects in 10-day global syntheses derived from VEGETATION/SPOT: II. Validation of an operational method on actual data sets [J]. Remote Sensing of Environment, 2002, 81(1): 101-113.

[4]       Huete, A. R., Didan, K., Miura, T., et al. Overview of the radiometric and biophysical performance of the MODIS vegetation indices [J]. Remote Sensing of Environment, 2002, 83(1/2): 195-213.

[5]       Leeuwen, W. J. D. V., Huete, A. R., Laing, T. W. MODIS vegetation index compositing approach: a prototype with AVHRR data [J]. Remote Sensing of Environment, 1999, 69(99): 264-280.

[6]       Leeuwen, W. J. D. V., Huete, A. R., Jia, S., et al. Comparison of vegetation index compositing scenarios: BRDF versus maximum VI approaches [C]. IEEE Geoscience and Remote Sensing Symposium, 1996, 1423-1425.

[7]       Walthall, C. L., Norman, J. M., Welles, J. M., et al. Simple equation to approximate the bidirectional reflectance from vegetative canopies and bare soil surfaces [J]. Applied Optics, 1985, 24(3): 383-387.

[8]       Vancutsem, C., Pekel, J. F., Bogaert, P., et al. Mean compositing, an alternative strategy for producing temporal syntheses [J]. International Journal of Remote Sensing, 2007, 28(22): 5123-5141.